CS180 - Project 2: Fun with Filters and Frequencies!

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import skimage as sk
import skimage.io as skio
from scipy.signal import convolve2d
import cv2
import os

Part 1.1: Convolutions from Scratch!

Note: What can you use for this section? This section is meant to be done with numpy only, simple array operations.

First, let's recap what a convolution is. Implement it with four for loops, then two for loops. Implement padding, with zero fill values; convolution without padding will receive partial credit.

Convolution formula:

$G[i, j] = \sum_{u=-k}^k \sum_{v=-k}^k H[u, v]F[i-u, j-v]$,

where $F$ is the image, $H$ is the kernel of size $2k+1$, and $G$ is the output image.

$kernelSize = 2k + 1 \implies k = \frac{kernelSize - 1}{2}$.

def convolution_four_loops(F, H):
    '''Convolve matrix, F, with kernel, H.'''
    k_x = (H.shape[0] - 1)//2
    k_y = (H.shape[1] - 1)//2
    F_padded = np.pad(F, ((k_x, k_x), (k_y, k_y)), mode='constant', constant_values=0)
    G = np.empty(F.shape)
    for i in range(0, G.shape[0]):
        for j in range(0, G.shape[1]):
            sum = 0
            for u in range(-k_x, k_x+1):
                for v in range(-k_y, k_y+1):
                    sum += H[u+k_x, v+k_y] * F_padded[i - u + k_x, j - v + k_y]
            G[i, j] = sum
    
    return G

def convolution_two_loops(F, H):
    '''Convolve matrix, F, with kernel, H.'''
    H_flipped = np.flip(H, axis=(0, 1)) #flip along horizontal and vertical axis to match four loops convolution
    k_x = (H.shape[0] - 1)//2
    k_y = (H.shape[1] - 1)//2
    F_padded = np.pad(F, ((k_x, k_x), (k_y, k_y)), mode='constant', constant_values=0)
    G = np.empty(F.shape)
    for i in range(0, G.shape[0]):
        for j in range(0, G.shape[1]):
            sum = 0
            region = F_padded[i:i+H.shape[0], j:j+H.shape[1]]
            sum = np.sum(region * H_flipped)
            G[i, j] = sum
    
    return G

Testing convolution implementation with finite difference operators $D_x$ and $D_y$

A = np.array([[0, 0, 0, 0, 0],
     [0, 0, 1, 0, 0],
     [0, 1, 1, 1, 0],
     [0, 0, 1, 0, 0],
     [0, 0, 0, 0, 0]])
D_x = np.array([[1, 0, -1]])
D_y = np.array([[1], [0], [-1]])

convolve_four_loops_x = convolution_four_loops(A, D_x)
convolve_four_loops_y = convolution_four_loops(A, D_y)

convolve_two_loops_x = convolution_two_loops(A, D_x)
convolve_two_loops_y = convolution_two_loops(A, D_y)

plt.figure(figsize=(10, 8))

# First row: only the original image
plt.subplot(3, 2, 1)
plt.imshow(A, cmap='gray')
plt.title('Original Image')
plt.axis('off')

# Second row: two images
plt.subplot(3, 2, 3)
plt.imshow(convolve_four_loops_x, cmap='gray')
plt.title('Convolved with D_x (4 loops)')
plt.axis('off')

plt.subplot(3, 2, 4)
plt.imshow(convolve_four_loops_y, cmap='gray')
plt.title('Convolved with D_y (4 loops)')
plt.axis('off')

# Third row: two images
plt.subplot(3, 2, 5)
plt.imshow(convolve_two_loops_x, cmap='gray')
plt.title('Convolved with D_x (2 loops)')
plt.axis('off')

plt.subplot(3, 2, 6)
plt.imshow(convolve_two_loops_y, cmap='gray')
plt.title('Convolved with D_y (2 loops)')
plt.axis('off')

plt.tight_layout()
plt.show()

Compare it with a built-in convolution function scipy.signal.convolve2d! Then, take a picture of yourself (and read it as grayscale), write out a 9x9 box filter, and convolve the picture with the box filter. Do it with the finite difference operators Dx and Dy as well. Include the code snippets in the website!

selfie = cv2.imread("data\joel_pic.jpg", cv2.IMREAD_GRAYSCALE)

# center = 1300
# width = 500
cropped_selfie = selfie #[center:center+width+1, center:center+width+1]

box_filter = np.ones((9,9)) / (81)

box_filtered_selfie_four_loops = convolution_four_loops(cropped_selfie, box_filter)

plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.imshow(cropped_selfie, cmap='gray')
plt.title('Selfie (Grayscale)')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(box_filtered_selfie_four_loops, cmap='gray')
plt.title('Box Filtered (Manual; four loops)') 
plt.axis('off')

plt.show()

box_filtered_selfie_two_loops = convolution_two_loops(cropped_selfie, box_filter)

plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.imshow(cropped_selfie, cmap='gray')
plt.title('Selfie (Grayscale)')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(box_filtered_selfie_two_loops, cmap='gray')
plt.title('Box Filtered (Manual; two loops)') 
plt.axis('off')

plt.show()

box_filtered_selfie_scipy = convolve2d(cropped_selfie, box_filter, mode='same', boundary='fill', fillvalue=0)

plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.imshow(cropped_selfie, cmap='gray')
plt.title('Selfie (Grayscale)')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(box_filtered_selfie_scipy, cmap='gray')
plt.title('Box Filtered (Scipy)') 
plt.axis('off')

plt.show()

In my two manual implementations of convolution, boundaries were handled in the same way: the input image was padded with 0s around the borders until it achieved the dimensions needed for a convolution with a kernel to be the same dimensions as the input image. This achieves the same effect as the mode I set the SciPy implementation to (ie. oundary='fill', fillvalue=0). It took about 12 minutes with my four loops manual approach, ~1 minute with two loops, and Scypy cracked it about 7 seconds.

The outputs looks pretty similar to the original. Next, lets take the difference of the filtered images with the original to see what was actually effected.

difference_manual_four_loops = cropped_selfie - box_filtered_selfie_four_loops
difference_manual_two_loops = cropped_selfie - box_filtered_selfie_two_loops
difference_scypy = cropped_selfie - box_filtered_selfie_scipy

plt.figure(figsize=(10, 5))
plt.subplot(1, 4, 1)
plt.imshow(cropped_selfie, cmap='gray')
plt.title('Selfie (Grayscale)')
plt.axis('off')

plt.subplot(1, 4, 2)
plt.imshow(difference_manual_four_loops, cmap='gray')
plt.title('Difference (four loops)')
plt.axis('off')

plt.subplot(1, 4, 3)
plt.imshow(difference_manual_two_loops, cmap='gray')
plt.title('Difference (two loops)')
plt.axis('off')

plt.subplot(1, 4, 4)
plt.imshow(difference_scypy, cmap='gray')
plt.title('Difference (scypy)')
plt.axis('off')
plt.tight_layout()

plt.show()

Upon closer inspection it is clear that the box filtered removed high frequency changes, most notably details in the hair. Zoom in and you'll notice its a bit more blurry accross that region. Further, all 3 implementations, despite the runtime differences, appear to be doing the achieving the same result (correctness looks promising; runtime could use some work!).

Next, let's convolve the selfie with the finite difference operators $D_x$ and $D_y$.

selfie_dx = convolution_two_loops(cropped_selfie, D_x)
selfie_dy = convolution_two_loops(cropped_selfie, D_y)
plt.figure(figsize=(12, 5))
plt.subplot(1, 3, 1)
plt.imshow(cropped_selfie, cmap='gray')
plt.title('Original (Grayscale)')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(selfie_dx, cmap='gray')
plt.title('Partial Derivative (Dx)')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(selfie_dx, cmap='gray')
plt.title('Partial Derivative (Dy)')
plt.axis('off')

plt.show()

Part 1.2: Finite Difference Operator

cameraman = cv2.imread("data\cameraman.png", cv2.IMREAD_GRAYSCALE)

First, show the partial derivative in x and y of the cameraman image by convolving the image with finite difference operators D_x and D_y.

cameraman_dx = convolve2d(cameraman, D_x, mode='same', boundary='fill', fillvalue=0)
cameraman_dy = convolve2d(cameraman, D_y, mode='same', boundary='fill', fillvalue=0)

plt.figure(figsize=(12, 5))
plt.subplot(1, 3, 1)
plt.imshow(cameraman, cmap='gray')
plt.title('Original (Grayscale)')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(cameraman_dx, cmap='gray')
plt.title('Partial Derivative (Dx)')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(cameraman_dy, cmap='gray')
plt.title('Partial Derivative (Dy)')
plt.axis('off')

plt.show()

Now compute and show the gradient magnitude image.

cameraman_gradient_mag = np.sqrt(cameraman_dx**2 + cameraman_dy**2)

plt.figure(figsize=(6, 6))
plt.imshow(cameraman_gradient_mag, cmap='gray')
plt.title('Gradient Magnitude')
plt.axis('off')
plt.show()

To turn this into an edge image, lets binarize the gradient magnitude image by picking the appropriate threshold (trying to suppress the noise while showing all the real edges; it will take you a few tries to find the right threshold; This threshold is meant to be assessed qualitatively).

cameraman_img_width, cameraman_img_height = cameraman_gradient_mag.shape
binarized_cameraman = np.zeros(cameraman_gradient_mag.shape)
threshold = 100
for x in range(cameraman_img_width):
    for y in range(cameraman_img_height):
        if cameraman_gradient_mag[x][y] > threshold:
            binarized_cameraman[x][y] = 255
        else:
            binarized_cameraman[x][y] = 0

plt.figure(figsize=(6, 6))
plt.imshow(binarized_cameraman, cmap='gray')
plt.title('Binarized Gradient Magnitude (pixel intensity threshold = 100)')
plt.axis('off')
plt.show()

The threshold for which values were high enough to be cieled (and which are floored) was set to 100 out 255. At this value, there is minimal "salt" around the image (i.e stray white pixels), yet the image is outline of the man and the camera is still primairly in tact. Just a bit higher a threshold and too much was floored (black image), and a bit lower results in a lot of noise making it through. A threshold of 100 was the sweat spot.

Part 1.3: Derivative of Gaussian (DoG) Filter

We noted that the results with just the difference operator were rather noisy. Luckily, we have a smoothing operator handy: the Gaussian filter G. Create a blurred version of the original image by convolving with a gaussian and repeat the procedure in the previous part (one way to create a 2D gaussian filter is by using cv2.getGaussianKernel() to create a 1D gaussian and then taking an outer product with its transpose to get a 2D gaussian kernel).

gaussian = cv2.getGaussianKernel(7, 1)
gaussian = gaussian * gaussian.T
cameraman_blurred = convolve2d(cameraman, gaussian, mode='same', boundary='fill', fillvalue=0)

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.imshow(cameraman, cmap='gray')
plt.title('Original (Grayscale)')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(cameraman_blurred, cmap='gray')
plt.title('Guassian Blur')
plt.axis('off')

(np.float64(-0.5), np.float64(539.5), np.float64(541.5), np.float64(-0.5))

cameraman_blurred_dx = convolve2d(cameraman_blurred, D_x, mode='same', boundary='fill', fillvalue=0)
cameraman_blurred_dy = convolve2d(cameraman_blurred, D_y, mode='same', boundary='fill', fillvalue=0)

cameraman_blurred_gradient_mag = np.sqrt(cameraman_blurred_dx**2 + cameraman_blurred_dy**2)

cameraman_blurred_img_width, cameraman_blurred_img_height = cameraman_blurred_gradient_mag.shape
binarized_cameraman_blurred_t100 = np.zeros(cameraman_blurred_gradient_mag.shape)
binarized_cameraman_blurred_t60 = np.zeros(cameraman_blurred_gradient_mag.shape)
for x in range(cameraman_blurred_img_width):
    for y in range(cameraman_blurred_img_height):
        if cameraman_blurred_gradient_mag[x][y] > 100: #threshold 1
            binarized_cameraman_blurred_t100[x][y] = 255
        else:
            binarized_cameraman_blurred_t100[x][y] = 0

        if cameraman_blurred_gradient_mag[x][y] > 60: #threshold 2
            binarized_cameraman_blurred_t60[x][y] = 255
        else:
            binarized_cameraman_blurred_t60[x][y] = 0

plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(binarized_cameraman_blurred_t100, cmap='gray')
plt.title('Gradient Mag (from blurred image, T=100)')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(binarized_cameraman_blurred_t60, cmap='gray')
plt.title('Gradient Mag (from blurred image, T=60)')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(binarized_cameraman, cmap='gray')
plt.title('Gradient Mag (from NOT blurred, T=100)')
plt.axis('off')

plt.subplots_adjust(wspace=0.5)

plt.show()

Check out the difference! I had to drop the threshold that we used to binarize the non-blurred gradient maginitude from 100 to 60 since the average intensity has been brought down; but, as you can see, the gradient magnitude from the blurred image at its best, produces much smoother boundary lines then the version produced from the non-blurred image.

Now we can do the same thing with a single convolution instead of two by creating a derivative of gaussian filters. Convolve the gaussian with D_x and D_y and display the resulting DoG filters as images.

Then, lets verify that we get the same result as before.

DoGx = convolve2d(gaussian, D_x, mode='same', boundary='fill', fillvalue=0)
DoGy = convolve2d(gaussian, D_y, mode='same', boundary='fill', fillvalue=0)

plt.figure(figsize=(10, 5))

# Plot DoGx (Derivative in X direction)
plt.subplot(1, 2, 1)
# Use 'seismic' or 'coolwarm' colormap to clearly show positive (red) and negative (blue) values
plt.imshow(DoGx, cmap='seismic')
plt.title(r'DoGx Kernel ($\frac{\partial G}{\partial x}$)')
plt.axis('off')
plt.colorbar(fraction=0.046, pad=0.04) # Add colorbar to show value range

# Plot DoGy (Derivative in Y direction)
plt.subplot(1, 2, 2)
plt.imshow(DoGy, cmap='seismic')
plt.title(r'DoGy Kernel ($\frac{\partial G}{\partial y}$)')
plt.axis('off')
plt.colorbar(fraction=0.046, pad=0.04) # Add colorbar

plt.tight_layout()
plt.show()

cameraman_DoGx = convolve2d(cameraman, DoGx, mode='same', boundary='fill', fillvalue=0)
cameraman_DoGy = convolve2d(cameraman, DoGy, mode='same', boundary='fill', fillvalue=0)

cameraman_DoG_gradient_mag = np.sqrt(cameraman_DoGx**2 + cameraman_DoGy**2)

cameraman_blurred_img_width, cameraman_blurred_img_height = cameraman_DoG_gradient_mag.shape
binarized_cameraman_DoG_t60 = np.zeros(cameraman_DoG_gradient_mag.shape)
for x in range(cameraman_blurred_img_width):
    for y in range(cameraman_blurred_img_height):
        if cameraman_DoG_gradient_mag[x][y] > 60:
            binarized_cameraman_DoG_t60[x][y] = 255
        else:
            binarized_cameraman_DoG_t60[x][y] = 0

plt.figure(figsize=(15, 5))

plt.subplot(1, 2, 1)
plt.imshow(binarized_cameraman_DoG_t60, cmap='gray')
plt.title('Gradient Mag (using DoG, T=60)')
plt.axis('off')


plt.subplot(1, 2, 2)
plt.imshow(binarized_cameraman_blurred_t60, cmap='gray')
plt.title('Gradient Mag (from blurred image, T=60)')
plt.axis('off')

plt.subplots_adjust(wspace=0.5)

plt.show()

And the resulting images are about the same (give or take a few pixels due to rounding errors between the two order in which to perform the image filtering), as desired.

Part 2: Fun with Frequencies!

Part 2.1: Image "Sharpening"

Pick your favorite blurry image and get ready to "sharpen" it! We will derive the unsharp masking technique. Remember our favorite Gaussian filter from class. This is a low pass filter that retains only the low frequencies. We can subtract the blurred version from the original image to get the high frequencies of the image. An image often looks sharper if it has stronger high frequencies. So, lets add a little bit more high frequencies to the image! Combine this into a single convolution operation which is called the unsharp mask filter.

1. Extract the low frequencies: $I_{\text{low}} = I \ast G$, where $I$ is the input image and $G$ is the gaussian kernal we used before to blur.
2. Extract the high frequencies: $I_{\text{high}} = I - I_{\text{low}}$
3. Add more of these high frequencies back to the original image: $I_{\text{sharpened}} = I + I_{\text{high}}$
   Optional: Lets add a scalar multiplier to $I_{\text{high}}$ to control how much more or less of these high frequencies we add: $I_{\text{sharpened}} = I + \alpha I_{\text{high}}$

To reduce these operations down to a single convoltion, lets rollout $I_{\text{sharpened}}$ to see how we can reduce the expression:
$I_{\text{sharpened}} = I + \alpha [I - I_{\text{low}}] = I + \alpha [I - I \ast G]$.
Note that $M \ast \delta = M$ for some matrix $M$, where $\delta$ is a square kernel that's all 0s, except for a 1 at its center.
$\implies I_{\text{sharpened}} = I \ast \delta + \alpha [ I \ast \delta - I \ast G]$
$\implies I_{\text{sharpened}} = I \ast \delta + I \ast \alpha[\delta - G]$, by commutativity of convolution
$\implies I_{\text{sharpened}} = I \ast [\delta + \alpha(\delta - G)]$.
Let $H_{\text{unsharp}} = \delta + \alpha(\delta - G)$.
$\implies I_{\text{sharpened}} = I \ast H_{\text{unsharp}}$, which is a single convolution as desired.
NOTE: I'm setting $\alpha = 10$ by default 'cause its is imperically a nice value to really notice the sharpening effects.

def get_Hunsharp(alpha=10):
    delta_kernel = np.zeros(gaussian.shape)
    delta_kernel[gaussian.shape[0]//2, gaussian.shape[1]//2] = 1.0

    H_unsharp = delta_kernel + alpha*(delta_kernel - gaussian)
    return H_unsharp

H_unsharp = get_Hunsharp(10)

Now to test on some images...

blurry_gray_cats = cv2.imread(r"data\blurry_cats.jpg", cv2.IMREAD_COLOR)
blurry_offended_cat = cv2.imread(r"data\offended_cat.jpg", cv2.IMREAD_COLOR)
blurry_taj_cat = cv2.imread(r"data\taj.jpg", cv2.IMREAD_COLOR)

blurry_gray_cats = cv2.cvtColor(blurry_gray_cats, cv2.COLOR_BGR2RGB)
blurry_offended_cat = cv2.cvtColor(blurry_offended_cat, cv2.COLOR_BGR2RGB)
blurry_taj_cat = cv2.cvtColor(blurry_taj_cat, cv2.COLOR_BGR2RGB)

if blurry_gray_cats is None or blurry_offended_cat is None or blurry_taj_cat is None:
    raise FileNotFoundError("One or more images failed to load. Check file paths and existence.")

def convolve_color_imgs(color_image, kernel, normalize = False, boundary='fill'):
    if color_image.dtype == np.uint8:
        I = color_image.astype(np.float32) / 255.0
    else:
        I = color_image.astype(np.float32)

    H, W, C = I.shape
    output_image = np.empty_like(I)

    for channel in range(C):
        single_channel = I[:, :, channel]

        if boundary == 'fill':
            filtered_channel = convolve2d(
                single_channel, 
                kernel, 
                mode='same', 
                boundary='fill', 
                fillvalue=0
            )
        elif boundary == 'symm':
            filtered_channel = convolve2d(
                single_channel,
                kernel,
                mode='same',
                boundary='symm'
            )

        output_image[:, :, channel] = filtered_channel

    if normalize:
        min_val = np.min(output_image)
        max_val = np.max(output_image)

        normalized_output_image = (output_image - min_val) / (max_val - min_val)
        output_image = normalized_output_image

    return output_image

sharp_gray_cats = convolve_color_imgs(blurry_gray_cats, H_unsharp)
sharp_offended_cat = convolve_color_imgs(blurry_offended_cat, H_unsharp)
sharp_taj_cat = convolve_color_imgs(blurry_taj_cat, H_unsharp)

plt.figure(figsize=(10, 15)) # Increased height for 3 rows

# --- ROW 1: Blurry Cats ---

# Left (Original/Blurry)
plt.subplot(3, 2, 1)
plt.imshow(blurry_gray_cats) 
plt.title('1. Blurry Cats (Original)')
plt.axis('off')

# Right (Sharpened)
plt.subplot(3, 2, 2)
plt.imshow(sharp_gray_cats)
plt.title('1. Blurry Cats (Sharpened)')
plt.axis('off')


# --- ROW 2: Offended Cat ---

# Left (Original/Blurry)
plt.subplot(3, 2, 3)
plt.imshow(blurry_offended_cat) 
plt.title('2. Offended Cat (Original)')
plt.axis('off')

# Right (Sharpened)
plt.subplot(3, 2, 4)
plt.imshow(sharp_offended_cat)
plt.title('2. Offended Cat (Sharpened)')
plt.axis('off')


# --- ROW 3: Taj Cat ---

# Left (Original/Blurry)
plt.subplot(3, 2, 5)
plt.imshow(blurry_taj_cat) 
plt.title('3. Taj Cat (Original)')
plt.axis('off')

# Right (Sharpened)
plt.subplot(3, 2, 6)
plt.imshow(sharp_taj_cat)
plt.title('3. Taj Cat (Sharpened)')
plt.axis('off')


# Increase vertical and horizontal spacing to prevent overlap
plt.subplots_adjust(wspace=0.1, hspace=0.15) 
plt.tight_layout()

plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.77370316..3.8300188].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.30125466..5.811322].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-4.8367715..4.834122].

Also for evaluation, lets pick a sharp image, blur it and then try to sharpen it again. Lets compare the original and the sharpened image:

blurred_sharp_taj_cat = convolve_color_imgs(sharp_taj_cat, gaussian)
sharp_blurred_sharp_taj_cat = convolve_color_imgs(blurred_sharp_taj_cat, H_unsharp)

plt.figure(figsize=(15, 5))

plt.subplot(1, 2, 1)
plt.imshow(blurry_taj_cat,)
plt.title('Taj cat (original)')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(sharp_blurred_sharp_taj_cat)
plt.title('Taj cat (original->sharpended->blurred->sharpened)')
plt.axis('off')

plt.subplots_adjust(wspace=0.5)

plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-9.169025..9.953941].

Now we're looking extra deep fried!

Finally, lets analyize how how the sharpening effect varies with our sharpness amount parameter $\alpha$:

alpha_values = list(range(1, 20, 3))
sharp_taj_cat_w_varied_alphas = []

for alpha in alpha_values:
    H_unsharp_alpha_tuned = get_Hunsharp(alpha)
    sharp_taj_cat_w_varied_alphas.append(convolve_color_imgs(blurry_taj_cat, H_unsharp_alpha_tuned))

fig, axes = plt.subplots(2, 4, figsize=(16, 12)) 
axes = axes.flatten() # Flatten the 3x4 array of axes for easy indexing

# Plot the Original Blurry Image first
axes[0].imshow(blurry_taj_cat)
axes[0].set_title(r'$\alpha=0$ (Original)', fontsize=12)
axes[0].axis('off')

# Iterate through the 11 sharpened images (starting at plot index 1)
for i, (alpha, sharp_img_raw) in enumerate(zip(alpha_values, sharp_taj_cat_w_varied_alphas)):
    # Calculate the plot index (start at 0 for the original, then 1, 2, 3...)
    plot_index = i + 1 
    
    # Normalize the calculated image to the [0, 1] range for proper display
    #sharpened_img_display = normalize_image_for_display(sharp_img_raw)
    
    axes[plot_index].imshow(sharp_img_raw)
    
    # Set the title using LaTeX for alpha (lambda)
    axes[plot_index].set_title(r'$\alpha={}$'.format(alpha), fontsize=12)
    axes[plot_index].axis('off')

# If the loop only produced 11 images (0 to 30 by 3), the last spot is empty.
# We can turn off the axis for the unused spot (index 11)
if len(sharp_taj_cat_w_varied_alphas) == 11:
    axes[11].axis('off')

plt.suptitle(r'Effect of Unsharp Mask $\alpha$ (Sharpening Amount) on Taj Cat Image', fontsize=16)
plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # Adjust for suptitle space
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.44838303..1.3396895].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-1.9111792..2.40659].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-3.3739753..3.620356].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-4.8367715..4.834122].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-6.2995677..6.0478883].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-7.7623634..7.2616544].
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-9.22516..8.47542].

Part 2.2: Hybrid Images

The goal of this part of the project is to create hybrid images using the approach described in the SIGGRAPH 2006 paper by Oliva, Torralba, and Schyns. Hybrid images are static images that change in interpretation as a function of the viewing distance. The basic idea is that high frequency tends to dominate perception when it is available, but, at a distance, only the low frequency (smooth) part of the signal can be seen. By blending the high frequency portion of one image with the low-frequency portion of another, you get a hybrid image that leads to different interpretations at different distances

1. First, we'll get a few pairs of images to make into hybrid images. Then, we will need to write code to low-pass filter one image, high-pass filter the second image, and add (or average) the two images. For a low-pass filter, Oliva et al. suggest using a standard 2D Gaussian filter. For a high-pass filter, they suggest using the impulse filter minus the Gaussian filter (which can be computed by subtracting the Gaussian-filtered image from the original). The cutoff-frequency of each filter should be chosen with some experimentation.

def get_delta_kernal_like(matrix):
    H, W = matrix.shape
    delta_kernel = np.zeros((H, W), dtype=np.float32)
    center_h, center_w = H // 2, W // 2
    
    # Set the center element to 1.0
    delta_kernel[center_h, center_w] = 1.0
    return delta_kernel

def make_gaussian_kernel(sigma):
    ksize = int(4 * sigma + 1)
    gaussian = cv2.getGaussianKernel(ksize, sigma)
    return gaussian @ gaussian.T

def color_hybrid_image(im1, im2, sigma1, sigma2, boundary='fill'):
    im1_low_pass = convolve_color_imgs(im1, make_gaussian_kernel(sigma1), boundary=boundary) # convolve with gaussian will result in just the low frequencies
    im2_low_pass = convolve_color_imgs(im2, make_gaussian_kernel(sigma2), boundary=boundary) # subtracting convolution with causing will leave just the higher frequencies
    im2_high_pass = im2 - im2_low_pass
    hybrid_raw = im1_low_pass + im2_high_pass
    hybrid_normalized = np.clip(hybrid_raw, 0, 1)
    return hybrid_normalized

def hybrid_image(im1, im2, sigma1, sigma2):
    if len(im1.shape) == 3:
        im1 = cv2.cvtColor(im1.astype(np.float32), cv2.COLOR_BGR2GRAY)
    if len(im2.shape) == 3:
        im2 = cv2.cvtColor(im2.astype(np.float32), cv2.COLOR_BGR2GRAY)

    # Normalize to [0,1]
    if im1.max() > 1.0:  
        im1 = im1 / 255.0
    if im2.max() > 1.0:
        im2 = im2 / 255.0

    gaussian1 = cv2.getGaussianKernel(7, sigma1)
    gaussian1 = gaussian1 @ gaussian1.T
    gaussian2 = cv2.getGaussianKernel(7, sigma2)
    gaussian2 = gaussian2 @ gaussian2.T

    im1_low_pass = convolve2d(im1, gaussian1, mode='same', boundary='fill', fillvalue=0)
    im2_low_pass = convolve2d(im2, gaussian2, mode='same', boundary='fill', fillvalue=0)
    im2_high_pass = im2 - im2_low_pass

    hybrid_raw = im1_low_pass + im2_high_pass
    
    hybrid_normalized = np.clip(hybrid_raw, 0, 1)
    return hybrid_normalized

%matplotlib qt
from align_image_code import align_images

# First load images
# high sf
derek = plt.imread(r'data/DerekPicture.jpg')/255.
# low sf
nutmeg = plt.imread(r'data/nutmeg.jpg')/255

# Next align images (this code is provided, but may be improved)
derek_aligned, nutmeg_aligned = align_images(derek, nutmeg)

## You will provide the code below. Sigma1 and sigma2 are arbitrary 
## cutoff values for the high and low frequencies
sigma1, sigma2 = 15, 15
hybrid = color_hybrid_image(derek_aligned, nutmeg_aligned, sigma1, sigma2)

%matplotlib inline

Please select 2 points in each image for alignment.

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

# First row: Derek and Nutmeg
axes[0, 0].imshow(derek)
axes[0, 0].set_title("Derek")
axes[0, 0].axis("off")

axes[0, 1].imshow(nutmeg)
axes[0, 1].set_title("Nutmeg")
axes[0, 1].axis("off")

# Second row: Hybrid image
axes[1, 0].imshow(hybrid)
axes[1, 0].set_title(f'Hybrid: Derek (Low Pass, $\\sigma_1={sigma1}$) + Nutmeg (High Pass, $\\sigma_2={sigma2}$)')
axes[1, 0].axis("off")

# Hide the empty subplot (bottom-right corner)
axes[1, 1].axis("off")

plt.tight_layout()
plt.show()

2. Finally, lets create 2-3 hybrid images (change of expression, morph between different objects, change over time, etc.). We'll show the input image and hybrid result per example.

%matplotlib qt
from align_image_code import align_images

# First load images
# high sf
cat_wizard = plt.imread(r'data/cat_wizard.png')
# low sf
joel = plt.imread(r'data/joel_pic.jpg')/255

# Next align images (this code is provided, but may be improved)
joel_aligned, cat_wizard_aligned = align_images(joel[::4, ::4], cat_wizard[::4, ::4])

## You will provide the code below. Sigma1 and sigma2 are arbitrary 
## cutoff values for the high and low frequencies
sigma1, sigma2 = 15, 15
hybrid2 = color_hybrid_image(joel_aligned, cat_wizard_aligned, sigma1, sigma2)

%matplotlib inline

Please select 2 points in each image for alignment.

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

# First row: Derek and Nutmeg
axes[0, 0].imshow(joel)
axes[0, 0].set_title("Joel")
axes[0, 0].axis("off")

axes[0, 1].imshow(cat_wizard)
axes[0, 1].set_title("Cat Wizard")
axes[0, 1].axis("off")

# Second row: Hybrid image
axes[1, 0].imshow(hybrid2)
axes[1, 0].set_title(f'Hybrid: Joel (Low Pass, $\\sigma_1={sigma1}$) + Cat Wizard (High Pass, $\\sigma_2={sigma2}$)')
axes[1, 0].axis("off")

# Hide the empty subplot (bottom-right corner)
axes[1, 1].axis("off")

plt.tight_layout()
plt.show()

%matplotlib qt
from align_image_code import align_images

# First load images
# high sf
doofenshmirtz = plt.imread(r'data/dr_doofenshmirtz.png')
# low sf
redestro = plt.imread(r'data/redestro.png')

# Next align images (this code is provided, but may be improved)
redestro_aligned, doofenshmirtz_aligned = align_images(redestro, doofenshmirtz)

## You will provide the code below. Sigma1 and sigma2 are arbitrary 
## cutoff values for the high and low frequencies
sigma1, sigma2 = 8, 10
hybrid3 = color_hybrid_image(redestro_aligned, doofenshmirtz_aligned, sigma1, sigma2)

%matplotlib inline

Please select 2 points in each image for alignment.

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

# First row: Derek and Nutmeg
axes[0, 0].imshow(doofenshmirtz)
axes[0, 0].set_title("Doofenshmirtz")
axes[0, 0].axis("off")

axes[0, 1].imshow(redestro)
axes[0, 1].set_title("Redestro")
axes[0, 1].axis("off")

# Second row: Hybrid image
axes[1, 0].imshow(hybrid3)
axes[1, 0].set_title(f'Hybrid: Doofenshmirtz (Low Pass, $\\sigma_1={sigma1}$) + Redestro (High Pass, $\\sigma_2={sigma2}$)')
axes[1, 0].axis("off")

# Hide the empty subplot (bottom-right corner)
axes[1, 1].axis("off")

plt.tight_layout()
plt.show()

3. For our favorite result, lets also illustrate the process through frequency analysis. We'll show the log magnitude of the Fourier transform of the two input images, the filtered images, and the hybrid image. We can use Python, to compute and display the 2D Fourier transform with: plt.imshow(np.log(np.abs(np.fft.fftshift(np.fft.fft2(gray_image)))))

def show_fft(image, title):
    gray = np.mean(image, axis=2) if image.ndim == 3 else image
    fft = np.fft.fftshift(np.fft.fft2(gray))
    magnitude = np.log(np.abs(fft) + 1e-6)  # add epsilon to avoid log(0)
    plt.imshow(magnitude, cmap="gray")
    plt.title(title)
    plt.axis("off")

im1_low_pass = convolve_color_imgs(redestro_aligned, make_gaussian_kernel(sigma1)) # convolve with gaussian will result in just the low frequencies
im2_low_pass = convolve_color_imgs(doofenshmirtz_aligned, make_gaussian_kernel(sigma2)) # subtracting convolution with causing will leave just the higher frequencies
im2_high_pass = doofenshmirtz_aligned - im2_low_pass
hybrid_raw = im1_low_pass + im2_high_pass
hybrid_normalized = np.clip(hybrid_raw, 0, 1)

fig, axes = plt.subplots(3, 4, figsize=(16, 12))

# Row 1: Original images + spectra
axes[0, 0].imshow(doofenshmirtz)
axes[0, 0].set_title("Doofenshmirtz")
axes[0, 0].axis("off")

plt.sca(axes[0, 1])
show_fft(doofenshmirtz, "FFT Doofenshmirtz")

axes[0, 2].imshow(redestro)
axes[0, 2].set_title("Redestro")
axes[0, 2].axis("off")

plt.sca(axes[0, 3])
show_fft(redestro, "FFT Redestro")

# Row 2: Filtered versions
axes[1, 0].imshow(im1_low_pass)
axes[1, 0].set_title("Red Low-pass")
axes[1, 0].axis("off")

plt.sca(axes[1, 1])
show_fft(im1_low_pass, "FFT Red Low-pass")

axes[1, 2].imshow(im2_high_pass + 0.5)  # shift for visibility
axes[1, 2].set_title("Doof High-pass")
axes[1, 2].axis("off")

plt.sca(axes[1, 3])
show_fft(im2_high_pass, "FFT Doof High-pass")

# Row 3: Hybrid result
axes[2, 0].imshow(hybrid3)
axes[2, 0].set_title("Hybrid (Doof + Red)")
axes[2, 0].axis("off")

plt.sca(axes[2, 1])
show_fft(hybrid3, "FFT Hybrid")

# Hide last two subplots in row 3
axes[2, 2].axis("off")
axes[2, 3].axis("off")

plt.tight_layout()
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.39986545..1.339786].

Part 2.3: Gaussian and Laplacian Stacks

In this part we will implement Gaussian and Laplacian stacks, which are kind of like pyramids but without the downsampling. This will prepare you for the next step for Multi-resolution blending.

1. Implement a Gaussian and a Laplacian stack. The different between a stack and a pyramid is that in each level of the pyramid the image is downsampled, so that the result gets smaller and smaller. In a stack the images are never downsampled so the results are all the same dimension as the original image, and can all be saved in one 3D matrix (if the original image was a grayscale image). To create the successive levels of the Gaussian Stack, just apply the Gaussian filter at each level, but do not subsample. In this way we will get a stack that behaves similarly to a pyramid that was downsampled to half its size at each level. If you would rather work with pyramids, you may implement pyramids other than stacks. However, in any case, you are NOT allowed to use built-in pyramid functions like cv2.pyrDown() or skimage.transform.pyramid_gaussian() in this project. You must implement your stacks from scratch!

#Modified from project 1 pyramid implementation
class ImageStack:
    def __init__(self, image, num_level=5, sigma=1.0):
        self.image = image
        self.num_levels = num_level
        self.sigma = sigma
        self.gaussian_stack = []
        self.laplacian_stack = []
        self._build_stacks()

    def _build_stacks(self):
        #first gaussian stack
        for level in range(self.num_levels):
            gaussian = make_gaussian_kernel(self.sigma * (level + 1))
            blurred_img = convolve_color_imgs(self.image, gaussian, boundary='symm')
            self.gaussian_stack.append(blurred_img)

        #then laplacian stack
        for i in range(self.num_levels - 1):
            laplacian_img = self.gaussian_stack[i] - self.gaussian_stack[i+1]
            self.laplacian_stack.append(laplacian_img)

        #set end of laplacian stack to the the highest frequency layer (most low frequencies removed) (i.e. the last level in gassian stack)
        self.laplacian_stack.append(self.gaussian_stack[-1])

    def get_gaussian_stack(self):
        return self.gaussian_stack

    def get_laplacian_stack(self):
        return self.laplacian_stack

apple = plt.imread(r'data/apple.png')
orange = plt.imread(r'data/orange.png')

apple = plt.imread(r'data/apple.png')
orange = plt.imread(r'data/orange.png')

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

ax[0].imshow(apple)
ax[0].set_title("Apple")
ax[0].axis("off")

ax[1].imshow(orange)
ax[1].set_title("Orange")
ax[1].axis("off")

plt.tight_layout()
plt.show()

apple_stack = ImageStack(apple, num_level=5, sigma=5.0)
orange_stack = ImageStack(orange, num_level=5, sigma=5.0)

apple_gaussians = apple_stack.get_gaussian_stack()
orange_gaussians = orange_stack.get_gaussian_stack()

fig, axes = plt.subplots(2, 5, figsize=(15, 6))

for i in range(5):
    axes[0, i].imshow(apple_gaussians[i])
    axes[0, i].set_title(f"Apple G{i}")
    axes[0, i].axis("off")

for i in range(5):
    axes[1, i].imshow(orange_gaussians[i])
    axes[1, i].set_title(f"Orange G{i}")
    axes[1, i].axis("off")

plt.tight_layout()
plt.show()

apple_laps = apple_stack.get_laplacian_stack()
orange_laps = orange_stack.get_laplacian_stack()

fig, axes = plt.subplots(2, 5, figsize=(15, 6))

for i in range(5):
    lap = apple_laps[i]
    lap_vis = (lap - lap.min()) / (lap.max() - lap.min())
    apple_laps[i] = lap_vis
    axes[0, i].imshow(lap_vis)
    axes[0, i].set_title(f"Apple L{i}")
    axes[0, i].axis("off")

for i in range(5):
    lap = orange_laps[i]
    lap_vis = (lap - lap.min()) / (lap.max() - lap.min())
    orange_laps[i] = lap_vis
    axes[1, i].imshow(lap_vis)
    axes[1, i].set_title(f"Orange L{i}")
    axes[1, i].axis("off")

plt.tight_layout()
plt.show()

2. Apply our Gaussian and Laplacian stacks to the Oraple and recreate the outcomes of Figure 3.42 in Szelski (Ed 2) page 167, as you can see in the image above. Review the 1983 paper for more information.

h, w, c = apple.shape
mask = np.zeros((h, w, c), dtype=np.float32)
mask[:, : w // 2, :] = 1.0  # left half 1 => apple
mask_stack = ImageStack(mask, num_level=5, sigma=15.0)

mask_gaussians = mask_stack.get_gaussian_stack()

fig, axes = plt.subplots(1, 5, figsize=(15, 6))

for i in range(5):
    axes[i].imshow(mask_gaussians[i])
    axes[i].set_title(f"Mask G{i}")
    axes[i].axis("off")

plt.tight_layout()
plt.show()

mask_gaussians = mask_stack.get_gaussian_stack()
num_levels = len(mask_gaussians)

fig, axes = plt.subplots(num_levels, 3, figsize=(12, 3*num_levels))

L_blend = []
for i, gaussian in enumerate(mask_gaussians):
    L_apple = apple_laps[i]
    L_orange = orange_laps[i]

    # ensure mask broadcasts if grayscale
    if gaussian.ndim == 2 and L_apple.ndim == 3:
        gaussian = gaussian[..., None]

    apple_contrib = gaussian * L_apple
    orange_contrib = (1 - gaussian) * L_orange
    L = apple_contrib + orange_contrib  # same as what we stored in L_blend
    L_blend.append(L)

    # plot
    axes[i, 0].imshow(np.clip(apple_contrib, 0, 1))
    axes[i, 0].set_title(f"Apple contrib L{i}")
    axes[i, 0].axis("off")

    axes[i, 1].imshow(np.clip(orange_contrib, 0, 1))
    axes[i, 1].set_title(f"Orange contrib L{i}")
    axes[i, 1].axis("off")

    axes[i, 2].imshow(np.clip(L, 0, 1))
    axes[i, 2].set_title(f"Blended L{i}")
    axes[i, 2].axis("off")

plt.tight_layout()
plt.show()

Part 2.4: Multiresolution Blending (a.k.a. the oraple!)

1. First, we'll need to get a few pairs of images that you want blend together with a vertical or horizontal seam. Then we will need to write some code in order to use our Gaussian and Laplacian stacks from part 2 in order to blend the images together. Since we are using stacks instead of pyramids like in the paper, the algorithm described on page 226 will not work as-is. If we try it out, you will find that you end up with a very clear seam between the apple and the orange since in the pyramid case the downsampling/blurring/upsampling hoopla ends up blurring the abrupt seam proposed in this algorithm. Instead, we should always use a mask as is proposed in the algorithm on page 230, and remember to create a Gaussian stack for your mask image as well as for the two input images. The Gaussian blurring of the mask in the pyramid will smooth out the transition between the two images. For the vertical or horizontal seam, our mask will simply be a step function of the same size as the original images.

Oraple = np.sum(L_blend, axis=0) # sum Laplacian levels
Oraple = np.clip(Oraple, 0, 1)

plt.figure(figsize=(6,6))
plt.imshow(Oraple)
plt.axis("off")
plt.title("Blended Oraple")
plt.show()

#weaker_H_unsharp = get_Hunsharp(5)
Oraple_sharp = convolve_color_imgs(convolve_color_imgs(Oraple, H_unsharp), H_unsharp)
plt.figure(figsize=(6,6))
plt.imshow(Oraple_sharp)
plt.axis("off")
plt.title("Blended Oraple")
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-2.6413076..28.855553].

Sharpened the image because it looked like some of our higher frequencies were diminished along the way. Looks great!

2. Now that we've made yourself an oraple (a.k.a your vertical or horizontal seam is nicely working), lets pick two pairs of images to blend together with an irregular mask, as is demonstrated in figure 8 in the paper.

andor_prison = plt.imread(r'data/andor_prison.png')
robot = plt.imread(r'data/joel_w_robot.png')
robot_mask = plt.imread(r'data/joel_w_robot_mask.png')

fig, ax = plt.subplots(1, 3, figsize=(10, 5))

ax[0].imshow(andor_prison)
ax[0].set_title("Andor Prison")
ax[0].axis("off")

ax[1].imshow(robot)
ax[1].set_title("Joel + Robot")
ax[1].axis("off")

ax[2].imshow(robot_mask)
ax[2].set_title("Joel + Robot (Mask)")
ax[2].axis("off")

plt.tight_layout()
plt.show()

andor_prison_stack = ImageStack(andor_prison, num_level=5, sigma=5.0)
robot_stack = ImageStack(robot, num_level=5, sigma=5.0)

andor_gaussians = andor_prison_stack.get_gaussian_stack()
robot_gaussians = robot_stack.get_gaussian_stack()

andor_lap = andor_prison_stack.get_laplacian_stack()
robot_lap = robot_stack.get_laplacian_stack()

for i in range(5):
    axes[0, i].imshow(andor_gaussians[i])
    axes[0, i].set_title(f"Andor Prison G{i}")
    axes[0, i].axis("off")

for i in range(5):
    axes[1, i].imshow(robot_gaussians[i])
    axes[1, i].set_title(f"Robot + Joel G{i}")
    axes[1, i].axis("off")

plt.tight_layout()
plt.show()

robot_mask_stack = ImageStack(robot_mask, num_level=5, sigma=5.0)

mask_gaussians = robot_mask_stack.get_gaussian_stack()
num_levels = len(mask_gaussians)
fig, axes = plt.subplots(1, num_levels, figsize=(3*num_levels, 3))

for i, g in enumerate(mask_gaussians):
    # If mask is single-channel (H,W), keep it grayscale
    if g.ndim == 2:
        axes[i].imshow(g, cmap="gray")
    else:  # if mask somehow has color channels
        axes[i].imshow(np.clip(g, 0, 1))
    axes[i].set_title(f"Mask G{i}")
    axes[i].axis("off")

plt.tight_layout()
plt.show()

mask_gaussians = robot_mask_stack.get_gaussian_stack()
num_levels = len(mask_gaussians)

Andor_L_blend = []
for i, gaussian in enumerate(mask_gaussians):
    L_left = andor_lap[i]
    L_right = robot_lap[i]

    # ensure mask broadcasts if grayscale
    if gaussian.ndim == 2 and L_apple.ndim == 3:
        gaussian = gaussian[..., None]

    L = gaussian * L_right + (1 - gaussian) * L_left  # same as what we stored in L_blend
    Andor_L_blend.append(L)

Andor_L_blend.append((robot_mask_stack.image * robot_stack.image + (1 - robot_mask_stack.image) * andor_prison_stack.image)*0.3)

andor_joel = np.sum(Andor_L_blend, axis=0) # sum Laplacian levels
andor_joel = np.clip(andor_joel, 0, 1)

plt.figure(figsize=(6,6))
plt.imshow(andor_joel)
plt.axis("off")
plt.title("Andor Prison + Joel Blended")
plt.show()

3. Finally, lets illustrate the process by applying your Laplacian stack and displaying it for your favorite result and the masked input images that created it. This should look similar to Figure 10 in the paper.

tito = plt.imread(r'data/tito.png')
sky = plt.imread(r'data/sky.png')
tito_mask = plt.imread(r'data/tito_mask.png')

fig, ax = plt.subplots(1, 3, figsize=(10, 5))

ax[0].imshow(tito)
ax[0].set_title("Andor Prison")
ax[0].axis("off")

ax[1].imshow(sky)
ax[1].set_title("Joel + Robot")
ax[1].axis("off")

ax[2].imshow(tito_mask)
ax[2].set_title("Joel + Robot (Mask)")
ax[2].axis("off")

plt.tight_layout()
plt.show()

tito_stack = ImageStack(tito, num_level=5, sigma=5.0)
sky_stack = ImageStack(sky, num_level=5, sigma=5.0)
tito_mask_stack = ImageStack(tito_mask, num_level=5, sigma=5.0)

tito_gaussians = tito_stack.get_gaussian_stack()
sky_gaussians = sky_stack.get_gaussian_stack()
tito_mask_gaussians = tito_mask_stack.get_gaussian_stack()

tito_lap = tito_stack.get_laplacian_stack()
sky_lap = sky_stack.get_laplacian_stack()
tito_mask_lap = tito_mask_stack.get_laplacian_stack()

fig, axes = plt.subplots(3, 5, figsize=(15, 9))  # 3 rows, 5 columns

for i in range(5):
    axes[0, i].imshow(tito_gaussians[i])
    axes[0, i].set_title(f"Tito G{i}")
    axes[0, i].axis("off")

for i in range(5):
    axes[1, i].imshow(sky_gaussians[i])
    axes[1, i].set_title(f"Sky G{i}")
    axes[1, i].axis("off")

for i in range(5):
    axes[2, i].imshow(tito_mask_gaussians[i])
    axes[2, i].set_title(f"Tito Mask G{i}")
    axes[2, i].axis("off")

plt.tight_layout()
plt.show()

mask_gaussians = tito_mask_gaussians
num_levels = len(mask_gaussians)

tito_L_blend = []
for i, gaussian in enumerate(mask_gaussians):
    L_left = tito_lap[i]
    L_right = sky_lap[i]

    # ensure mask broadcasts if grayscale
    if gaussian.ndim == 2 and L_apple.ndim == 3:
        gaussian = gaussian[..., None]

    L = gaussian * L_right + (1 - gaussian) * L_left  # same as what we stored in L_blend
    tito_L_blend.append(L)

tito_L_blend.append((tito_mask_stack.image * sky_stack.image + (1 - tito_mask_stack.image)* tito_stack.image)*0.3)

tito_sky = np.sum(tito_L_blend, axis=0) # sum Laplacian levels
tito_sky = np.clip(tito_sky, 0, 1)

plt.figure(figsize=(6,6))
plt.imshow(tito_sky)
plt.axis("off")
plt.title("Tito + Sky")
plt.show()

mask_gaussians = tito_mask_gaussians
num_levels = len(mask_gaussians)

tito_L_blend = []
for i, gaussian in enumerate(mask_gaussians):
    L_left = sky_lap[i]
    L_right = tito_lap[i]

    # ensure mask broadcasts if grayscale
    if gaussian.ndim == 2 and L_apple.ndim == 3:
        gaussian = gaussian[..., None]

    L = gaussian * L_right + (1 - gaussian) * L_left  # same as what we stored in L_blend
    tito_L_blend.append(L)

tito_L_blend.append((tito_mask_stack.image * tito_stack.image+ (1 - tito_mask_stack.image)* sky_stack.image)*0.3)

tito_sky = np.sum(tito_L_blend, axis=0) # sum Laplacian levels
tito_sky = np.clip(tito_sky, 0, 1)

plt.figure(figsize=(6,6))
plt.imshow(tito_sky)
plt.axis("off")
plt.title("Sky + Tito")
plt.show()